home *** CD-ROM | disk | FTP | other *** search
/ Ham Radio 2000 #2 / Ham Radio 2000 - Volume 2.iso / HAMV2 / MISC / COIL200 / MATHL.C < prev    next >
C/C++ Source or Header  |  1996-04-16  |  7KB  |  369 lines
  1. #include <stdlib.h>
  2. #include <stdio.h>
  3.  
  4. extern double MACHEP, MAXNUM;
  5.  
  6. /* Factorials of integers from 0 through 33 */
  7. double fac[] = {
  8.   1.0,
  9.   1.0,
  10.   2.0,
  11.   6.0,
  12.   2.4E1,
  13.   1.2E2,
  14.   7.2E2,
  15.   5.04E3,
  16.   4.032E4,
  17.   3.6288E5,
  18.   3.6288E6,
  19.   3.99168E7,
  20.   4.790016E8,
  21.   6.2270208E9,
  22.   8.71782912E10,
  23.   1.307674368E12,
  24.   2.0922789888E13,
  25.   3.55687428096E14,
  26.   6.402373705728E15,
  27.   1.21645100408832E17,
  28.   2.43290200817664E18,
  29.   5.109094217170944E19,
  30.   1.12400072777760768E21,
  31.   2.585201673888497664E22,
  32.   6.2044840173323943936E23,
  33.   1.5511210043330985984E25,
  34.   4.03291461126605635584E26,
  35.   1.0888869450418352160768E28,
  36.   3.04888344611713860501504E29,
  37.   8.841761993739701954543616E30,
  38.   2.6525285981219105863630848E32,
  39.   8.22283865417792281772556288E33,
  40.   2.6313083693369353016721801216E35,
  41.   8.68331761881188649551819440128E36
  42. };
  43.  
  44.  
  45. /*                            ellpe.c
  46.  *
  47.  *    Complete elliptic integral of the second kind
  48.  *
  49.  *
  50.  *
  51.  * SYNOPSIS:
  52.  *
  53.  * double m1, y, ellpe();
  54.  *
  55.  * y = ellpe( m1 );
  56.  *
  57.  *
  58.  *
  59.  * DESCRIPTION:
  60.  *
  61.  * Approximates the integral
  62.  *
  63.  *
  64.  *          pi/2
  65.  *           -
  66.  *          | |          2
  67.  * E(m)     =    |       sqrt( 1 - m sin t ) dt
  68.  *        | |
  69.  *         -
  70.  *          0
  71.  *
  72.  * Where m = 1 - m1, using the approximation
  73.  *
  74.  *    P(x)  -     x log x Q(x).
  75.  *
  76.  * Though there are no singularities, the argument m1 is used
  77.  * rather than m for compatibility with ellpk().
  78.  *
  79.  * E(1) = 1; E(0) = pi/2.
  80.  *
  81.  *
  82.  * ACCURACY:
  83.  *
  84.  *            Relative error:
  85.  * arithmetic    domain       # trials     peak          rms
  86.  *    DEC     0, 1        13000    3.1e-17        9.4e-18
  87.  *    IEEE     0, 1        10000    2.1e-16        7.3e-17
  88.  *
  89.  *
  90.  * ERROR MESSAGES:
  91.  *
  92.  *   message         condition        value returned
  93.  * ellpe domain         x<0, x>1         0.0
  94.  *
  95.  */
  96.  
  97. /*
  98. Cephes Math Library, Release 2.1:  February, 1989
  99. Copyright 1984, 1987, 1989 by Stephen L. Moshier
  100. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  101. */
  102.  
  103. static double PE[] = {
  104.   1.53552577301013293365E-4,
  105.   2.50888492163602060990E-3,
  106.   8.68786816565889628429E-3,
  107.   1.07350949056076193403E-2,
  108.   7.77395492516787092951E-3,
  109.   7.58395289413514708519E-3,
  110.   1.15688436810574127319E-2,
  111.   2.18317996015557253103E-2,
  112.   5.68051945617860553470E-2,
  113.   4.43147180560990850618E-1,
  114.   1.00000000000000000299E0
  115. };
  116. static double QE[] = {
  117.   3.27954898576485872656E-5,
  118.   1.00962792679356715133E-3,
  119.   6.50609489976927491433E-3,
  120.   1.68862163993311317300E-2,
  121.   2.61769742454493659583E-2,
  122.   3.34833904888224918614E-2,
  123.   4.27180926518931511717E-2,
  124.   5.85936634471101055642E-2,
  125.   9.37499997197644278445E-2,
  126.   2.49999999999888314361E-1
  127. };
  128.  
  129. double ellpe(x)
  130. double x;
  131. {
  132. double polevl(), log();
  133.  
  134.  
  135. if( (x <= 0.0) || (x > 1.0) )
  136.     {
  137.     if( x == 0.0 )
  138.         return( 1.0 );
  139.     printf( "ellpe domain error\n" );
  140.     return( 0.0 );
  141.     }
  142.  
  143. return( polevl(x,PE,10) - log(x) * (x * polevl(x,QE,9)) );
  144. }
  145.  
  146. /*                            ellpk.c
  147.  *
  148.  *    Complete elliptic integral of the first kind
  149.  *
  150.  *
  151.  *
  152.  * SYNOPSIS:
  153.  *
  154.  * double m1, y, ellpk();
  155.  *
  156.  * y = ellpk( m1 );
  157.  *
  158.  *
  159.  *
  160.  * DESCRIPTION:
  161.  *
  162.  * Approximates the integral
  163.  *
  164.  *
  165.  *
  166.  *          pi/2
  167.  *           -
  168.  *          | |
  169.  *          |          dt
  170.  * K(m)     =    |       ------------------
  171.  *          |              2
  172.  *        | |       sqrt( 1 - m sin t )
  173.  *         -
  174.  *          0
  175.  *
  176.  * where m = 1 - m1, using the approximation
  177.  *
  178.  *     P(x)  -    log x Q(x).
  179.  *
  180.  * The argument m1 is used rather than m so that the logarithmic
  181.  * singularity at m = 1 will be shifted to the origin; this
  182.  * preserves maximum accuracy.
  183.  *
  184.  * K(0) = pi/2.
  185.  *
  186.  * ACCURACY:
  187.  *
  188.  *            Relative error:
  189.  * arithmetic    domain       # trials     peak          rms
  190.  *    DEC     0,1        16000    3.5e-17        1.1e-17
  191.  *    IEEE     0,1        30000    2.5e-16        6.8e-17
  192.  *
  193.  * ERROR MESSAGES:
  194.  *
  195.  *   message         condition        value returned
  196.  * ellpk domain          x<0, x>1         0.0
  197.  *
  198.  */
  199.  
  200.  
  201. /*
  202. Cephes Math Library, Release 2.0:  April, 1987
  203. Copyright 1984, 1987 by Stephen L. Moshier
  204. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  205. */
  206.  
  207. static double PK[] =
  208. {
  209.  1.37982864606273237150E-4,
  210.  2.28025724005875567385E-3,
  211.  7.97404013220415179367E-3,
  212.  9.85821379021226008714E-3,
  213.  6.87489687449949877925E-3,
  214.  6.18901033637687613229E-3,
  215.  8.79078273952743772254E-3,
  216.  1.49380448916805252718E-2,
  217.  3.08851465246711995998E-2,
  218.  9.65735902811690126535E-2,
  219.  1.38629436111989062502E0
  220. };
  221.  
  222. static double QK[] =
  223. {
  224.  2.94078955048598507511E-5,
  225.  9.14184723865917226571E-4,
  226.  5.94058303753167793257E-3,
  227.  1.54850516649762399335E-2,
  228.  2.39089602715924892727E-2,
  229.  3.01204715227604046988E-2,
  230.  3.73774314173823228969E-2,
  231.  4.88280347570998239232E-2,
  232.  7.03124996963957469739E-2,
  233.  1.24999999999870820058E-1,
  234.  4.99999999999999999821E-1
  235. };
  236. static double C1 = 1.3862943611198906188E0; /* log(4) */
  237.  
  238.  
  239.  
  240. double ellpk(x)
  241. double x;
  242. {
  243. double polevl(), log();
  244.  
  245.  
  246.  
  247. if( (x < 0.0) || (x > 1.0) )
  248.     {
  249.     printf( "ellpk domain error\n" );
  250.     return( 0.0 );
  251.     }
  252.  
  253. if( x > MACHEP )
  254.     {
  255.     return( polevl(x,PK,10) - log(x) * polevl(x,QK,10) );
  256.     }
  257. else
  258.     {
  259.     if( x == 0.0 )
  260.         {
  261.         printf( "ellpk singularity\n" );
  262.         return( MAXNUM );
  263.         }
  264.     else
  265.         {
  266.         return( C1 - 0.5 * log(x) );
  267.         }
  268.     }
  269. }
  270.  
  271.  
  272. /*                            polevl.c
  273.  *                            p1evl.c
  274.  *
  275.  *    Evaluate polynomial
  276.  *
  277.  *
  278.  *
  279.  * SYNOPSIS:
  280.  *
  281.  * int N;
  282.  * double x, y, coef[N+1], polevl[];
  283.  *
  284.  * y = polevl( x, coef, N );
  285.  *
  286.  *
  287.  *
  288.  * DESCRIPTION:
  289.  *
  290.  * Evaluates polynomial of degree N:
  291.  *
  292.  *               2      N
  293.  * y  =     C  + C x + C x     +...+ C x
  294.  *      0    1     2        N
  295.  *
  296.  * Coefficients are stored in reverse order:
  297.  *
  298.  * coef[0] = C    , ..., coef[N] = C  .
  299.  *          N              0
  300.  *
  301.  *  The function p1evl() assumes that coef[N] = 1.0 and is
  302.  * omitted from the array.  Its calling arguments are
  303.  * otherwise the same as polevl().
  304.  *
  305.  *
  306.  * SPEED:
  307.  *
  308.  * In the interest of speed, there are no checks for out
  309.  * of bounds arithmetic.  This routine is used by most of
  310.  * the functions in the library.  Depending on available
  311.  * equipment features, the user may wish to rewrite the
  312.  * program in microcode or assembly language.
  313.  *
  314.  */
  315.  
  316.  
  317. /*
  318. Cephes Math Library Release 2.1:  December, 1988
  319. Copyright 1984, 1987, 1988 by Stephen L. Moshier
  320. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  321. */
  322.  
  323.  
  324. double polevl( x, coef, N )
  325. double x;
  326. double coef[];
  327. int N;
  328. {
  329. double ans;
  330. int i;
  331. double *p;
  332.  
  333. p = coef;
  334. ans = *p++;
  335. i = N;
  336.  
  337. do
  338.     ans = ans * x  +  *p++;
  339. while( --i );
  340.  
  341. return( ans );
  342. }
  343.  
  344. /*                            p1evl() */
  345. /*                        N
  346.  * Evaluate polynomial when coefficient of x  is 1.0.
  347.  * Otherwise same as polevl.
  348.  */
  349.  
  350. double p1evl( x, coef, N )
  351. double x;
  352. double coef[];
  353. int N;
  354. {
  355. double ans;
  356. double *p;
  357. int i;
  358.  
  359. p = coef;
  360. ans = x + *p++;
  361. i = N-1;
  362.  
  363. do
  364.     ans = ans * x  + *p++;
  365. while( --i );
  366.  
  367. return( ans );
  368. }
  369.